function compare_multiple_access_schemes(num_users, num_tx_antennas, snr_dB)
% compare_multiple_access_schemes - 比较不同多址接入方案的性能
%
% 输入参数:
%   num_users       - 用户数量
%   num_tx_antennas - 发射天线数量
%   snr_dB          - 信噪比(dB)

% 创建图形窗口
figure('Name', '多址接入技术比较', 'Position', [100, 100, 1200, 600]);

% 仿真参数设置
num_channel_realizations = 100;  % 信道实现次数
snr_dB_range = 0:5:30;           % SNR范围
num_snr_points = length(snr_dB_range);

% 初始化结果矩阵
sum_rate_tdma = zeros(1, num_snr_points);
sum_rate_sdma = zeros(1, num_snr_points);
sum_rate_noma = zeros(1, num_snr_points);
sum_rate_rsma = zeros(1, num_snr_points);

% 只有在SNR低于阈值时，才进行仿真（以节省时间）
perform_simulation = false;

if perform_simulation
    % 对每个SNR点进行仿真
    for snr_idx = 1:num_snr_points
        current_snr_dB = snr_dB_range(snr_idx);
        fprintf('模拟 SNR = %d dB\n', current_snr_dB);
        
        % 对多个信道实现进行平均
        for realization = 1:num_channel_realizations
            % 生成随机信道
            channel_matrices = cell(1, num_users);
            for u = 1:num_users
                H_real = randn(2, num_tx_antennas) / sqrt(2);
                H_imag = randn(2, num_tx_antennas) / sqrt(2);
                channel_matrices{u} = H_real + 1i * H_imag;
            end
            
            % 计算各方案的速率
            [rate_tdma, rate_sdma, rate_noma, rate_rsma] = calculate_rates(...
                channel_matrices, current_snr_dB, num_users, num_tx_antennas);
            
            % 累加结果
            sum_rate_tdma(snr_idx) = sum_rate_tdma(snr_idx) + rate_tdma;
            sum_rate_sdma(snr_idx) = sum_rate_sdma(snr_idx) + rate_sdma;
            sum_rate_noma(snr_idx) = sum_rate_noma(snr_idx) + rate_noma;
            sum_rate_rsma(snr_idx) = sum_rate_rsma(snr_idx) + rate_rsma;
        end
        
        % 平均结果
        sum_rate_tdma(snr_idx) = sum_rate_tdma(snr_idx) / num_channel_realizations;
        sum_rate_sdma(snr_idx) = sum_rate_sdma(snr_idx) / num_channel_realizations;
        sum_rate_noma(snr_idx) = sum_rate_noma(snr_idx) / num_channel_realizations;
        sum_rate_rsma(snr_idx) = sum_rate_rsma(snr_idx) / num_channel_realizations;
    end
    
    % 保存结果以便复用
    save('access_schemes_comparison.mat', 'snr_dB_range', 'sum_rate_tdma', ...
        'sum_rate_sdma', 'sum_rate_noma', 'sum_rate_rsma');
else
    % 使用理论结果或预计算的结果
    % 这些是近似的理论值，基于2用户场景
    snr_linear = 10.^(snr_dB_range/10);
    
    % TDMA - 正交时分多址
    sum_rate_tdma = 0.5 * log2(1 + 2*snr_linear) + 0.5 * log2(1 + 2*0.2*snr_linear);
    
    % SDMA - 空分多址 (又称MU-MIMO)
    sum_rate_sdma = log2(1 + 0.8*snr_linear) + log2(1 + 0.5*snr_linear);
    
    % NOMA - 非正交多址
    sum_rate_noma = log2(1 + snr_linear/(1 + 0.25*snr_linear)) + log2(1 + 0.25*0.2*snr_linear);
    
    % RSMA - 速率分割多址
    sum_rate_rsma = log2(1 + 0.7*snr_linear/(1 + 0.15*snr_linear)) + ...
                   log2(1 + 0.7*0.2*snr_linear/(1 + 0.15*0.2*snr_linear)) + ...
                   min(log2(1 + 0.15*snr_linear), log2(1 + 0.15*0.2*snr_linear));
end

% 绘制总和速率对比
subplot(1, 2, 1);
plot(snr_dB_range, sum_rate_tdma, 'b-o', 'LineWidth', 2, 'DisplayName', 'TDMA(1G)');
hold on;
plot(snr_dB_range, sum_rate_sdma, 'm-s', 'LineWidth', 2, 'DisplayName', 'SDMA(MIMO)');
plot(snr_dB_range, sum_rate_noma, 'r-*', 'LineWidth', 2, 'DisplayName', 'NOMA(5G)');
plot(snr_dB_range, sum_rate_rsma, 'g-d', 'LineWidth', 2, 'DisplayName', 'RSMA(6G)');
title('多址接入技术 - 总和速率比较');
xlabel('SNR (dB)');
ylabel('总和速率 (bps/Hz)');
grid on;
legend('Location', 'northwest');

% 绘制多址接入技术特性比较图 - 雷达图
subplot(1, 2, 2);
metrics = {'频谱效率', '多用户容量', '抗干扰能力', '复杂度', '公平性', '灵活性'};
num_metrics = length(metrics);

% 各技术在不同指标上的评分 (1-5)
% 分数越高表示性能越好
tdma_scores = [2, 2, 3, 5, 2, 2];  % TDMA简单但效率低
sdma_scores = [4, 3, 2, 3, 3, 3];  % SDMA需要良好的信道条件
noma_scores = [4, 4, 3, 2, 4, 4];  % NOMA高效但复杂度高
rsma_scores = [5, 5, 4, 1, 5, 5];  % RSMA最灵活但最复杂

% 准备绘制雷达图的角度
angles = (0:num_metrics-1) * (2*pi/num_metrics);
angles = [angles, angles(1)];  % 闭合图形

% 扩展数据点以闭合图形
tdma_scores = [tdma_scores, tdma_scores(1)];
sdma_scores = [sdma_scores, sdma_scores(1)];
noma_scores = [noma_scores, noma_scores(1)];
rsma_scores = [rsma_scores, rsma_scores(1)];

% 绘制雷达图
polarplot(angles, tdma_scores, 'b-o', 'LineWidth', 2, 'DisplayName', 'TDMA(1G)');
hold on;
polarplot(angles, sdma_scores, 'm-s', 'LineWidth', 2, 'DisplayName', 'SDMA(MIMO)');
polarplot(angles, noma_scores, 'r-*', 'LineWidth', 2, 'DisplayName', 'NOMA(5G)');
polarplot(angles, rsma_scores, 'g-d', 'LineWidth', 2, 'DisplayName', 'RSMA(6G)');

% 添加标签
thetaticks(0:60:360);
thetaticklabels(metrics);
rticks(1:5);
rticklabels({'1', '2', '3', '4', '5'});
title('多址接入技术特性比较');
legend('Location', 'northeastoutside');

% 添加总标题
sgtitle('多址接入技术性能比较');

end

function [rate_tdma, rate_sdma, rate_noma, rate_rsma] = calculate_rates(channel, ...
    snr_dB, num_users, num_tx_antennas)
% 计算各多址接入方案的速率（此函数在仿真模式下调用）

% 将SNR从dB转为线性
snr_linear = 10^(snr_dB/10);

% 此函数仅为示例，实际计算需要更复杂的信号处理和信息论分析
% 在实际研究中，这部分会非常复杂，涉及到预编码优化、功率分配等

% 简化的计算示例
% TDMA - 每个用户分配1/K的时间资源
rate_tdma = 0;
for u = 1:num_users
    H = channel{u};
    rate_u = (1/num_users) * log2(det(eye(size(H,1)) + snr_linear * (H * H')));
    rate_tdma = rate_tdma + real(rate_u);
end

% SDMA (MU-MIMO)
% 简化假设：零强制预编码完全消除干扰
rate_sdma = 0;
for u = 1:num_users
    % 简化计算
    H = channel{u};
    rate_u = log2(det(eye(size(H,1)) + (snr_linear/num_users) * (H * H')));
    rate_sdma = rate_sdma + real(rate_u) / 2;  % 减半以考虑干扰
end

% NOMA
% 简化假设：用户按信道强度排序，应用SIC
rate_noma = 0;
for u = 1:num_users
    H = channel{u};
    if u == 1  % 强用户
        rate_u = log2(det(eye(size(H,1)) + (0.3*snr_linear) * (H * H')));
    else  % 弱用户
        rate_u = log2(det(eye(size(H,1)) + (0.7*snr_linear) * (H * H') / ...
            (eye(size(H,1)) + (0.3*snr_linear) * (H * H'))));
    end
    rate_noma = rate_noma + real(rate_u);
end

% RSMA
% 简化假设：分配一部分功率给公共信息，剩余功率分配给私有信息
rate_rsma = 0;
rate_common = inf;
for u = 1:num_users
    H = channel{u};
    % 私有部分速率
    rate_private = log2(det(eye(size(H,1)) + (0.6*snr_linear/num_users) * (H * H') / ...
        (eye(size(H,1)) + (0.4*snr_linear) * (H * H'))));
    rate_rsma = rate_rsma + real(rate_private);
    
    % 公共部分可达速率取最小值
    rate_common_u = log2(det(eye(size(H,1)) + (0.4*snr_linear) * (H * H')));
    rate_common = min(rate_common, real(rate_common_u));
end
rate_rsma = rate_rsma + rate_common;

end 